options(scipen=100,digits=3)
library(cmdstanr)
options(mc.cores=parallel::detectCores())
options(cmdstanr_max_rows=1000)
library(gridExtra)

execute mcmc sampling

goMCMC=function(mdl,data,smp=500,wrm=100,th=1){
  mcmc=mdl$sample(
  data=data,
  seed=1,
  chains=4,
  iter_sampling=smp,
  iter_warmup=wrm,
  thin=th,
  refresh=1000
  )
  mcmc
}

see mcmc result and parameters

seeMCMC=function(mcmc,exc='',ch=T){ # exclude 'exc' parameters from seeing
  print(mcmc)
  prs=mcmc$metadata()$model_params[-1] # reject lp__
  for(pr in prs){
    if(grepl('^y',pr)) next # not show predictive value "y*" information
    if(exc!='' && grepl(paste0('^',exc),pr)) next
    drw=mcmc$draws(pr)
    if(ch){
      par(mfrow=c(4,1),mar=c(1,5,1,1))
      drw[,1,] |> plot(type='l',xlab='',ylab=pr)
      drw[,2,] |> plot(type='l',xlab='',ylab=pr)
      drw[,3,] |> plot(type='l',xlab='',ylab=pr)
      drw[,4,] |> plot(type='l',xlab='',ylab=pr)
      par(mar=c(3,5,3,3))
    }

    par(mfrow=c(1,2))
    drw |> hist(main=pr,xlab='')
    drw |> density() |> plot(main=pr)    
  }
}

stan basic

normal distribution

ex1-0.stan

normal distribution

data{
  int N;
  vector[N] y;
}
parameters{
  real m;
  real<lower=0> s;
}
model{
  y~normal(m,s);
}
generated quantities{
  real y1;
  y1=normal_rng(m,s);
}
mdl=cmdstan_model('./ex1-0.stan')
y=rnorm(10,2,1)
summary(y)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    1.24    1.76    2.37    2.41    2.74    3.89
par(mfrow=c(1,2))
hist(y,main='y')
density(y) |> plot(main='y')  

data=list(N=length(y),y=y)

mle=mdl$optimize(data=data)
## Initial log joint probability = -10.039 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       11      -2.72838   5.37318e-05   2.33335e-07           1           1       15    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.3 seconds.
mle
##  variable estimate
##      lp__    -2.73
##      m        2.41
##      s        0.80
##      y1       4.14
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.3 seconds.
mcmc
##  variable  mean median   sd  mad    q5   q95 rhat ess_bulk ess_tail
##      lp__ -3.98  -3.67 1.01 0.75 -5.96 -3.00 1.00      913      998
##      m     2.41   2.42 0.30 0.28  1.91  2.90 1.00     1049     1011
##      s     0.98   0.93 0.28 0.24  0.65  1.50 1.00     1332     1265
##      y1    2.37   2.37 1.03 1.00  0.69  4.02 1.00     1726     1581
mcmc$metadata()$stan_variables
## [1] "lp__" "m"    "s"    "y1"
mcmc$metadata()$model_params
## [1] "lp__" "m"    "s"    "y1"
seeMCMC(mcmc)
##  variable  mean median   sd  mad    q5   q95 rhat ess_bulk ess_tail
##      lp__ -3.98  -3.67 1.01 0.75 -5.96 -3.00 1.00      913      998
##      m     2.41   2.42 0.30 0.28  1.91  2.90 1.00     1049     1011
##      s     0.98   0.93 0.28 0.24  0.65  1.50 1.00     1332     1265
##      y1    2.37   2.37 1.03 1.00  0.69  4.02 1.00     1726     1581

poisson distribution

y=rpois(20,3)
summary(y)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##     0.0     2.0     3.5     3.3     5.0     8.0
par(mfrow=c(1,2))
hist(y,main='y')
density(y) |> plot(main='y')

data=list(N=length(y),y=y)

ex2-1.stan

in case fit poisson dist. to normal dist.

data{
  int N;
  vector[N] y;
}

parameters{
  real<lower=0> m;
  real<lower=0> s;
}

model{
  y~normal(m,s);
}

generated quantities{
  real y1;
  y1=normal_rng(m,s);
}
# when fitting poisson dist. sample to normal dist.
mdl=cmdstan_model('./ex2-1.stan')

mle=mdl$optimize(data=data)
## Initial log joint probability = -345.95 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       14      -22.8371   0.000167347    0.00112883      0.9391      0.9391       16    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##      lp__   -22.84
##      m        3.30
##      s        1.90
##      y1       0.12
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc)
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##      lp__ -22.07 -21.74 1.10 0.80 -24.20 -21.04 1.00      878      826
##      m      3.28   3.28 0.49 0.46   2.47   4.07 1.00     1179      857
##      s      2.08   2.03 0.36 0.34   1.59   2.73 1.00     1198     1125
##      y1     3.31   3.32 2.22 2.16  -0.36   6.97 1.00     1904     1972

y1=mcmc$draws('y1')
par(mfrow=c(1,2))
hist(y1,main='y1')
density(y1) |> plot(main='y1')     

ex2-2.stan

fit to poisson dist.

data{
  int N;
  array[N] int y;
}

parameters{
  real<lower=0> l;
}

model{
  y~poisson(l);
}

generated quantities{
  int y1;
  y1=poisson_rng(l);
}
mdl=cmdstan_model('./ex2-2.stan')

mle=mdl$optimize(data=data)
## Initial log joint probability = -8.21616 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##        6       12.7989   0.000110146   1.60607e-05           1           1        7    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##      lp__    12.80
##      l        3.30
##      y1       1.00
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc)
##  variable  mean median   sd  mad    q5   q95 rhat ess_bulk ess_tail
##      lp__ 13.49  13.77 0.70 0.31 12.02 14.00 1.00     1172     1478
##      l     3.37   3.35 0.41 0.41  2.72  4.09 1.01      881     1161
##      y1    3.36   3.00 1.87 1.48  1.00  7.00 1.00     1881     1892

y1=mcmc$draws('y1')
par(mfrow=c(1,2))
hist(y1,main='y1')
density(y1) |> plot(main='y1') 

difference test

ex3.stan

mean difference of 2 normal distributions

data{
  int N;
  vector[N] a;
  vector[N] b;
}
parameters{
  real ma;
  real<lower=0> sa;
  real mb;
  real<lower=0> sb;
}
model{
  a~normal(ma,sa);
  b~normal(mb,sb);
}
generated quantities{
  real d;
  d=mb-ma;
}
n=20
a=rnorm(n,10,1)
b=rnorm(n,11,2)
data=list(N=n,a=a,b=b)

mdl=cmdstan_model('./ex3.stan')

mle=mdl$optimize(data=data)
## Initial log joint probability = -954.202 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       25      -33.0393   8.26231e-05   0.000276317           1           1       45    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##      lp__   -33.04
##      ma       9.70
##      sa       0.99
##      mb      11.09
##      sb       1.93
##      d        1.39
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.3 seconds.
seeMCMC(mcmc,ch=F)
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##      lp__ -34.47 -34.09 1.52 1.32 -37.48 -32.73 1.00      800      945
##      ma     9.71   9.70 0.25 0.23   9.31  10.13 1.00     1536     1192
##      sa     1.09   1.06 0.20 0.18   0.82   1.45 1.00     2343     1507
##      mb    11.09  11.10 0.48 0.45  10.28  11.88 1.00      844      899
##      sb     2.13   2.07 0.39 0.36   1.61   2.85 1.00     2097     1467
##      d      1.38   1.38 0.52 0.51   0.53   2.24 1.00      932     1027

d=mcmc$draws('d')
d
## # A draws_array: 500 iterations, 4 chains, and 1 variables
## , , variable = d
## 
##          chain
## iteration    1   2   3    4
##         1 1.99 1.8 1.9 1.30
##         2 0.84 1.7 2.1 1.04
##         3 1.32 1.9 1.4 1.00
##         4 0.74 1.5 1.9 0.58
##         5 1.43 1.0 2.0 1.01
## 
## # ... with 495 more iterations
mean(d>0)
## [1] 0.993

normal regression

ex4-1.stan

single normal regression

data{
  int N;
  vector[N] x;
  vector[N] y;
}

parameters{
  real b0;
  real b1;
  real<lower=0> s;
}

model{
  y~normal(b0+b1*x,s);
}
# make prediction with R
n=20
x=runif(n,0,100)
y=rnorm(n,x*3+10,10)
par(mfrow=c(1,1))
plot(x,y)

data=list(N=n,x=x,y=y)

mdl=cmdstan_model('./ex4-1.stan')

mle=mdl$optimize(data=data)
## Initial log joint probability = -3.43654e+07 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
## Error evaluating model log probability: Non-finite gradient. 
## Error evaluating model log probability: Non-finite gradient. 
## Error evaluating model log probability: Non-finite gradient. 
## Error evaluating model log probability: Non-finite gradient. 
## Error evaluating model log probability: Non-finite gradient. 
## Error evaluating model log probability: Non-finite gradient. 
## Error evaluating model log probability: Non-finite gradient. 
## Error evaluating model log probability: Non-finite gradient. 
## Error evaluating model log probability: Non-finite gradient. 
##       40      -59.6172    0.00359792   0.000762108           1           1       98    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##      lp__   -59.62
##      b0      13.84
##      b1       2.92
##      s       11.95
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 finished in 0.2 seconds.
## Chain 4 finished in 0.2 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.3 seconds.
seeMCMC(mcmc,ch=F)
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##      lp__ -58.68 -58.32 1.28 1.03 -61.24 -57.29 1.01      474      554
##      b0    13.93  14.08 6.22 5.94   3.36  23.95 1.01      298      261
##      b1     2.91   2.91 0.10 0.10   2.76   3.08 1.01      389      464
##      s     13.43  13.07 2.33 2.07  10.28  17.62 1.01     1125      943

b0=mcmc$draws('b0')
b1=mcmc$draws('b1')
b=tibble(b0=c(b0),b1=c(b1)) |> sample_n(100)
x1=runif(nrow(b),min(x),max(x))
xy=tibble(x=x1,y=b$b0+b$b1*x1)
par(mfrow=c(1,1))
plot(xy)

ex4-2.stan

single normal regression, prediction from new explanatory variable x1

data{
  int N;
  vector[N] x;
  vector[N] y;
  int N1;
  vector[N1] x1;
}

parameters{
  real b0;
  real b1;
  real<lower=0> s;
}

model{
  y~normal(b0+b1*x,s);
}

generated quantities{
  vector[N1] m;
  vector[N1] y1;
  for(i in 1:N1){
    m[i]=b0+b1*x1[i];
    y1[i]=normal_rng(m[i],s);
  }
}
# make prediction with stan
n1=10
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcition

data=list(N=n,x=x,y=y,N1=n1,x1=x1)

mdl=cmdstan_model('./ex4-2.stan')

mle=mdl$optimize(data=data)
## Initial log joint probability = -1.21236e+06 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
## Error evaluating model log probability: Non-finite gradient. 
## Error evaluating model log probability: Non-finite gradient. 
## Exception: normal_lpdf: Scale parameter is 0, but must be positive! (in '/tmp/RtmpSo844c/model-6216618232e.stan', line 16, column 2 to column 22) 
## Error evaluating model log probability: Non-finite gradient. 
## Exception: normal_lpdf: Scale parameter is 0, but must be positive! (in '/tmp/RtmpSo844c/model-6216618232e.stan', line 16, column 2 to column 22) 
## Exception: normal_lpdf: Scale parameter is 0, but must be positive! (in '/tmp/RtmpSo844c/model-6216618232e.stan', line 16, column 2 to column 22) 
## Exception: normal_lpdf: Scale parameter is 0, but must be positive! (in '/tmp/RtmpSo844c/model-6216618232e.stan', line 16, column 2 to column 22) 
## Error evaluating model log probability: Non-finite gradient. 
## Error evaluating model log probability: Non-finite gradient. 
##       42      -59.6172    0.00334236    0.00415075      0.8495      0.8495      101    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##    lp__     -59.62
##    b0        13.84
##    b1         2.92
##    s         11.95
##    m[1]      25.57
##    m[2]      55.51
##    m[3]      85.45
##    m[4]     115.39
##    m[5]     145.33
##    m[6]     175.27
##    m[7]     205.21
##    m[8]     235.15
##    m[9]     265.09
##    m[10]    295.03
##    y1[1]     37.14
##    y1[2]     71.44
##    y1[3]     80.16
##    y1[4]    120.42
##    y1[5]    146.37
##    y1[6]    163.83
##    y1[7]    206.59
##    y1[8]    229.99
##    y1[9]    277.56
##    y1[10]   293.50
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 finished in 0.1 seconds.
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 finished in 0.2 seconds.
## Chain 3 finished in 0.2 seconds.
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 finished in 0.2 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.2 seconds.
## Total execution time: 0.3 seconds.
mcmc$metadata()$stan_variables
## [1] "lp__" "b0"   "b1"   "s"    "m"    "y1"
seeMCMC(mcmc,'m',ch=F)
##  variable   mean median    sd   mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -58.60 -58.30  1.16  0.99 -60.96 -57.31 1.01      490      727
##    b0      14.01  14.02  5.82  5.67   4.84  23.33 1.02      321      264
##    b1       2.91   2.91  0.10  0.09   2.75   3.07 1.01      413      447
##    s       13.49  13.13  2.35  2.19  10.27  17.95 1.00      943      913
##    m[1]    25.73  25.75  5.49  5.39  17.08  34.51 1.02      321      255
##    m[2]    55.64  55.72  4.69  4.69  48.30  63.23 1.02      325      255
##    m[3]    85.56  85.56  3.98  3.94  79.27  92.01 1.02      342      299
##    m[4]   115.47 115.43  3.40  3.34 110.06 121.00 1.01      394      466
##    m[5]   145.39 145.36  3.03  3.02 140.32 150.38 1.00      550      650
##    m[6]   175.30 175.31  2.96  2.97 170.43 180.12 1.00     1007     1342
##    m[7]   205.21 205.31  3.20  3.11 199.82 210.35 1.00     1941     1585
##    m[8]   235.13 235.16  3.69  3.57 228.84 241.04 1.00     2020     1402
##    m[9]   265.04 265.15  4.36  4.21 257.66 271.96 1.00     1459     1473
##    m[10]  294.96 295.07  5.12  4.89 286.37 303.22 1.00     1104     1407
##    y1[1]   25.94  26.12 14.82 14.52   1.99  50.38 1.00     1313     1800
##    y1[2]   55.56  55.86 14.09 13.12  31.65  77.20 1.00     1472     1496
##    y1[3]   85.53  85.29 14.12 13.67  62.72 108.89 1.00     1610     1802
##    y1[4]  115.78 115.60 14.51 13.55  91.61 139.51 1.00     1450     1525
##    y1[5]  145.68 145.54 13.95 13.33 122.81 168.64 1.00     2048     1702
##    y1[6]  175.48 175.22 13.89 13.49 153.83 199.25 1.00     1980     1759
##    y1[7]  204.72 204.31 14.05 13.57 183.15 227.86 1.00     2048     1740
##    y1[8]  235.19 235.08 14.37 14.03 211.79 258.63 1.00     2068     1819
##    y1[9]  264.87 264.95 14.42 13.83 241.73 288.04 1.00     2211     1912
##    y1[10] 294.99 294.93 14.49 13.89 271.06 318.63 1.00     1600     1972

m=mcmc$draws('m')
summary(m)
## # A tibble: 10 × 10
##    variable  mean median    sd   mad    q5   q95  rhat ess_bulk ess_tail
##    <chr>    <dbl>  <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>    <dbl>    <dbl>
##  1 m[1]      25.7   25.7  5.49  5.39  17.1  34.5  1.02     321.     255.
##  2 m[2]      55.6   55.7  4.69  4.69  48.3  63.2  1.02     326.     256.
##  3 m[3]      85.6   85.6  3.98  3.94  79.3  92.0  1.02     342.     299.
##  4 m[4]     115.   115.   3.40  3.34 110.  121.   1.01     394.     467.
##  5 m[5]     145.   145.   3.03  3.02 140.  150.   1.00     551.     650.
##  6 m[6]     175.   175.   2.96  2.97 170.  180.   1.00    1008.    1342.
##  7 m[7]     205.   205.   3.20  3.11 200.  210.   1.00    1942.    1586.
##  8 m[8]     235.   235.   3.69  3.57 229.  241.   1.00    2020.    1403.
##  9 m[9]     265.   265.   4.36  4.21 258.  272.   1.00    1459.    1473.
## 10 m[10]    295.   295.   5.12  4.89 286.  303.   1.00    1104.    1407.
smm=summary(m)

y1=mcmc$draws('y1')
summary(y1)
## # A tibble: 10 × 10
##    variable  mean median    sd   mad     q5   q95  rhat ess_bulk ess_tail
##    <chr>    <dbl>  <dbl> <dbl> <dbl>  <dbl> <dbl> <dbl>    <dbl>    <dbl>
##  1 y1[1]     25.9   26.1  14.8  14.5   1.99  50.4  1.00    1314.    1800.
##  2 y1[2]     55.6   55.9  14.1  13.1  31.6   77.2  1.00    1473.    1497.
##  3 y1[3]     85.5   85.3  14.1  13.7  62.7  109.   1.00    1610.    1802.
##  4 y1[4]    116.   116.   14.5  13.5  91.6  140.   1.00    1450.    1525.
##  5 y1[5]    146.   146.   14.0  13.3 123.   169.   1.00    2048.    1702.
##  6 y1[6]    175.   175.   13.9  13.5 154.   199.   1.00    1981.    1759.
##  7 y1[7]    205.   204.   14.0  13.6 183.   228.   1.00    2049.    1740.
##  8 y1[8]    235.   235.   14.4  14.0 212.   259.   1.00    2069.    1820.
##  9 y1[9]    265.   265.   14.4  13.8 242.   288.   1.00    2211.    1912.
## 10 y1[10]   295.   295.   14.5  13.9 271.   319.   1.00    1600.    1973.
smy=summary(y1)

xy=tibble(x=x1,m=smm$median,ml5=smm$q5,mu5=smm$q95,yl5=smy$q5,yu5=smy$q95)

par(mfrow=c(1,1))
xlim=c(min(x),max(x))
ylim=c(min(y),max(y))
plot(x,y,
     xlim=xlim,ylim=ylim, xlab='x',ylab='y',col='red')
par(new=T)
plot(xy$x,xy$m,type='l',
     xlim=xlim,ylim=ylim, xlab='',ylab='',col='black')
par(new=T)
plot(xy$x,xy$ml5,type='l',
     xlim=xlim,ylim=ylim, xlab='',ylab='',col='darkgray')
par(new=T)
plot(xy$x,xy$mu5,type='l',
     xlim=xlim,ylim=ylim, xlab='',ylab='',col='darkgray')
par(new=T)
plot(xy$x,xy$yl5,type='l',
     xlim=xlim,ylim=ylim, xlab='',ylab='',col='lightgray')
par(new=T)
plot(xy$x,xy$yu5,type='l',
     xlim=xlim,ylim=ylim, xlab='',ylab='',col='lightgray')

qplot(x,y,col=I('red'))+
  geom_line(aes(x=x,y=m),xy,col='black')+
  geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
  geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
  geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
  geom_line(aes(x=x,y=yu5),xy,col='lightgray')

ex4-3.stan

single normal regression, prediction from original explanatory variable x

data{
  int N;
  vector[N] x;
  vector[N] y;
}

parameters{
  real b0;
  real b1;
  real<lower=0> s;
}

model{
  y~normal(b0+b1*x,s);
}

generated quantities{
  vector[N] m;
  vector[N] y1;
  for(i in 1:N){
    m[i]=b0+b1*x[i];
    y1[i]=normal_rng(m[i],s);
  }
}
data=list(N=n,x=x,y=y)

mdl=cmdstan_model('./ex4-3.stan')

mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 finished in 0.1 seconds.
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 finished in 0.2 seconds.
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 finished in 0.2 seconds.
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 finished in 0.2 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.2 seconds.
## Total execution time: 0.3 seconds.
seeMCMC(mcmc,'m',ch=F)
##  variable   mean median    sd   mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -58.66 -58.31  1.28  1.04 -61.17 -57.30 1.01      404      876
##    b0      13.31  13.13  5.87  5.47   3.85  22.77 1.02      177      202
##    b1       2.93   2.93  0.10  0.09   2.77   3.09 1.02      242      362
##    s       13.46  13.15  2.47  2.30  10.19  18.12 1.00     1002     1041
##    m[1]    69.11  69.02  4.37  4.17  62.00  76.37 1.02      193      285
##    m[2]   258.00 257.93  4.31  3.85 251.23 265.34 1.00     1393     1335
##    m[3]    71.90  71.80  4.30  4.11  64.91  78.94 1.02      196      295
##    m[4]    77.86  77.79  4.16  4.01  71.08  84.65 1.02      202      355
##    m[5]    25.07  24.91  5.53  5.18  16.04  34.03 1.02      179      205
##    m[6]   295.27 295.23  5.27  4.85 286.97 304.29 1.01      940      956
##    m[7]    96.33  96.25  3.76  3.60  90.18 102.55 1.02      231      426
##    m[8]    32.73  32.60  5.32  4.98  23.95  41.44 1.02      179      213
##    m[9]   279.40 279.33  4.85  4.36 271.98 287.70 1.00     1086     1076
##    m[10]  212.77 212.74  3.39  3.19 207.24 218.39 1.00     2087     1384
##    m[11]  182.41 182.36  3.04  2.86 177.48 187.40 1.00     1266     1448
##    m[12]  235.76 235.69  3.81  3.52 229.67 242.09 1.00     1911     1432
##    m[13]  181.29 181.25  3.04  2.84 176.36 186.29 1.00     1232     1447
##    m[14]  161.08 161.03  2.99  2.74 156.24 166.03 1.00      762      986
##    m[15]  194.13 194.11  3.14  2.98 188.98 199.40 1.00     1629     1365
##    m[16]  232.51 232.46  3.74  3.48 226.47 238.79 1.00     1936     1445
##    m[17]   76.23  76.14  4.20  4.04  69.40  83.07 1.02      200      316
##    m[18]  269.00 268.93  4.58  4.12 261.81 276.86 1.00     1215     1107
##    m[19]  145.28 145.22  3.06  2.82 140.35 150.34 1.01      552      897
##    m[20]  275.90 275.82  4.75  4.27 268.58 284.04 1.00     1127     1076
##    y1[1]   68.77  68.82 14.69 13.86  44.11  92.56 1.01     1101     1665
##    y1[2]  258.30 258.08 14.92 14.37 233.46 283.22 1.00     1916     1702
##    y1[3]   72.28  72.16 13.99 13.10  49.03  95.75 1.00     1788     1892
##    y1[4]   77.88  77.90 14.16 13.48  54.65 101.28 1.00     1465     1583
##    y1[5]   25.19  25.32 14.63 14.56   1.64  48.97 1.00     1359     1574
##    y1[6]  295.41 294.86 15.13 14.32 271.77 319.90 1.00     1726     1609
##    y1[7]   96.58  96.60 14.19 13.65  72.99 119.14 1.00     1631     1845
##    y1[8]   32.74  32.72 14.63 14.33   7.70  56.44 1.01     1017     1555
##    y1[9]  279.61 279.56 14.28 14.07 256.50 303.27 1.00     2008     1842
##    y1[10] 212.26 212.05 14.04 13.92 189.96 235.49 1.00     2145     1825
##    y1[11] 182.24 182.44 14.08 13.22 159.52 205.44 1.00     1951     1721
##    y1[12] 235.67 235.66 14.13 13.52 212.14 257.93 1.00     1993     1847
##    y1[13] 181.64 181.63 13.89 13.51 159.82 204.80 1.00     1946     1798
##    y1[14] 161.13 161.55 14.00 13.70 138.42 183.56 1.00     1983     1871
##    y1[15] 194.16 194.13 13.82 13.37 171.67 217.11 1.00     2040     2014
##    y1[16] 232.64 232.76 14.59 14.18 208.35 256.20 1.00     2156     1921
##    y1[17]  76.37  76.68 14.41 13.98  53.13  99.92 1.00     1692     1883
##    y1[18] 269.23 268.76 14.36 13.53 246.17 293.34 1.00     2016     2002
##    y1[19] 144.85 145.02 13.99 12.37 120.95 167.95 1.00     1875     1605
##    y1[20] 276.05 276.39 14.22 13.29 252.34 298.78 1.00     1684     1976

y1=mcmc$draws('y1')
smy=summary(y1)

par(mfrow=c(1,1))
plot(y,smy$median,xlab='obs.',ylab='prd.')
abline(0,1)

qplot(y,smy$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1)

par(mfrow=c(1,2))
hist(y-smy$median,xlab='obs.-prd.',main='residual')
density(y-smy$median) |> plot(xlab='obs.-prd.',main='residual')

estimate correlation

ex4-41.stan

use covariance matrix and multinormal distribution

data{
  int N;
  array[N] vector[2] y;
}

parameters{
  vector[2] m;
  vector<lower=0>[2] s;
  real<lower=-1,upper=1> r;
}

transformed parameters{
  cov_matrix[2] cv;
  cv[1,1]=s[1]^2;
  cv[2,2]=s[2]^2;
  cv[1,2]=r*s[1]*s[2];
  cv[2,1]=r*s[1]*s[2];
}

model{
  y~multi_normal(m,cv);
}

ex4-42.stan

use covariance matrix and wishart distribution

data{
  int N;
  matrix[2,2] dp;
}

parameters{
  vector<lower=0>[2] s;
  real<lower=-1,upper=1> r;
}

transformed parameters{
  cov_matrix[2] cv;
  cv[1,1]=s[1]^2;
  cv[2,2]=s[2]^2;
  cv[1,2]=r*s[1]*s[2];
  cv[2,1]=r*s[1]*s[2];
}

model{
  dp~wishart(N-1,cv);
}

ex4-43.stan

use correlation matrix and wishart distribution

data{
  int N;
  int M;
  matrix[M,M] dp;
}

parameters{
  vector<lower=0>[M] s;
  corr_matrix[M] r;
}

transformed parameters{
  cov_matrix[M] cv;
  cv=quad_form_diag(r,s);
}

model{
  dp~wishart(N-1,cv);
}

ex4-44.stan use covariance matrix and multinormal distribution

data{
  int N;
  int K;
  array[N] vector[K] y;
}
parameters{
  vector[K] m;
  cov_matrix[K] cov;
}
model{
  y~multi_normal(m,cov);
}
library(MASS)
n=20
y=mvrnorm(n,c(0,0),matrix(c(1,0.7,0.7,1),2),2)
cov(y)
##      [,1] [,2]
## [1,] 1.10 0.59
## [2,] 0.59 1.33
cor(y)
##       [,1]  [,2]
## [1,] 1.000 0.489
## [2,] 0.489 1.000
plot(y)

#estimate covariance
mdl=cmdstan_model('./ex4-41.stan')

data=list(N=n,y=y)

mle=mdl$optimize(data=data)
## Initial log joint probability = -185.469 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       16      -20.0064   0.000311034    6.1786e-05           1           1       18    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##   lp__      -20.01
##   m[1]        0.32
##   m[2]        0.44
##   s[1]        1.02
##   s[2]        1.12
##   r           0.49
##   cv[1,1]     1.04
##   cv[2,1]     0.56
##   cv[1,2]     0.56
##   cv[2,2]     1.26
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc,ch=F)
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##   lp__    -23.38 -23.04 1.59 1.44 -26.29 -21.42 1.00      965     1280
##   m[1]      0.33   0.33 0.26 0.25  -0.10   0.75 1.00     1405     1124
##   m[2]      0.44   0.44 0.29 0.28  -0.02   0.92 1.00     1365     1128
##   s[1]      1.15   1.12 0.20 0.19   0.87   1.51 1.00     2131     1602
##   s[2]      1.25   1.23 0.22 0.20   0.96   1.65 1.00     1717     1771
##   r         0.45   0.47 0.18 0.18   0.11   0.72 1.00      772     1131
##   cv[1,1]   1.36   1.25 0.49 0.41   0.76   2.29 1.00     2131     1602
##   cv[2,1]   0.68   0.62 0.42 0.34   0.14   1.44 1.00      826     1101
##   cv[1,2]   0.68   0.62 0.42 0.34   0.14   1.44 1.00      826     1101
##   cv[2,2]   1.62   1.50 0.61 0.48   0.91   2.73 1.00     1717     1771

#estimate correlation,covariance
#deviation product sum matirix folows Wishart dist.
#deviation product sum = covariance * n-1 or n

mdl=cmdstan_model('./ex4-42.stan')

data=list(N=n,dp=cov(y)*(n-1))

mle=mdl$optimize(data=data)
## Initial log joint probability = -40.578 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       11      -19.9806   0.000169014   0.000224214      0.9947      0.9947       14    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##   lp__      -19.98
##   s[1]        1.05
##   s[2]        1.15
##   r           0.49
##   cv[1,1]     1.10
##   cv[2,1]     0.59
##   cv[1,2]     0.59
##   cv[2,2]     1.33
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc,ch=F)
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##   lp__    -22.30 -21.96 1.26 1.02 -24.81 -20.93 1.00      741      980
##   s[1]      1.15   1.12 0.21 0.20   0.87   1.52 1.00     1608     1423
##   s[2]      1.27   1.24 0.23 0.23   0.95   1.70 1.00     1383     1264
##   r         0.45   0.47 0.18 0.18   0.12   0.73 1.00      622      485
##   cv[1,1]   1.36   1.26 0.53 0.44   0.75   2.31 1.00     1608     1423
##   cv[2,1]   0.69   0.62 0.43 0.36   0.16   1.51 1.00      686      679
##   cv[1,2]   0.69   0.62 0.43 0.36   0.16   1.51 1.00      686      679
##   cv[2,2]   1.66   1.53 0.63 0.55   0.90   2.88 1.00     1383     1264

#estimate correlation,covariance
mdl=cmdstan_model('./ex4-43.stan')

m=2
data=list(N=n,M=m,dp=cov(y)*(n-1))

mle=mdl$optimize(data=data)
## Initial log joint probability = -1495.83 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       15      -19.9806   3.96584e-05   0.000146643           1           1       19    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##   lp__      -19.98
##   s[1]        1.05
##   s[2]        1.15
##   r[1,1]      1.00
##   r[2,1]      0.49
##   r[1,2]      0.49
##   r[2,2]      1.00
##   cv[1,1]     1.10
##   cv[2,1]     0.59
##   cv[1,2]     0.59
##   cv[2,2]     1.33
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc,ch=F)
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##   lp__    -21.56 -21.22 1.27 1.03 -24.06 -20.19 1.00      712      995
##   s[1]      1.14   1.12 0.21 0.20   0.86   1.52 1.00     1187     1030
##   s[2]      1.26   1.23 0.23 0.20   0.96   1.68 1.00     1155      914
##   r[1,1]    1.00   1.00 0.00 0.00   1.00   1.00   NA       NA       NA
##   r[2,1]    0.45   0.47 0.19 0.18   0.12   0.73 1.00      850      835
##   r[1,2]    0.45   0.47 0.19 0.18   0.12   0.73 1.00      850      835
##   r[2,2]    1.00   1.00 0.00 0.00   1.00   1.00   NA       NA       NA
##   cv[1,1]   1.35   1.25 0.52 0.44   0.74   2.32 1.00     1187     1030
##   cv[2,1]   0.69   0.62 0.45 0.35   0.14   1.49 1.00      767      602
##   cv[1,2]   0.69   0.62 0.45 0.35   0.14   1.49 1.00      767      602
##   cv[2,2]   1.64   1.51 0.63 0.48   0.92   2.82 1.00     1155      914

mdl=cmdstan_model('./ex4-44.stan')

k=2
data=list(N=n,K=k,y=y)

mle=mdl$optimize(data=data)
## Initial log joint probability = -16515.4 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
## Exception: multi_normal_lpdf: LDLT_Factor of covariance parameter is not positive definite.  last conditional variance is -3.42114e-49. (in '/tmp/RtmpbEIVDP/model-5029291c0ba9.stan', line 11, column 2 to column 24) 
## Exception: multi_normal_lpdf: LDLT_Factor of covariance parameter is not positive definite.  last conditional variance is 0. (in '/tmp/RtmpbEIVDP/model-5029291c0ba9.stan', line 11, column 2 to column 24) 
## Exception: multi_normal_lpdf: LDLT_Factor of covariance parameter is not positive definite.  last conditional variance is -1.72723e-77. (in '/tmp/RtmpbEIVDP/model-5029291c0ba9.stan', line 11, column 2 to column 24) 
## Exception: multi_normal_lpdf: LDLT_Factor of covariance parameter is not positive definite.  last conditional variance is 0. (in '/tmp/RtmpbEIVDP/model-5029291c0ba9.stan', line 11, column 2 to column 24) 
## Exception: multi_normal_lpdf: LDLT_Factor of covariance parameter is not positive definite.  last conditional variance is 0. (in '/tmp/RtmpbEIVDP/model-5029291c0ba9.stan', line 11, column 2 to column 24) 
##       44      -20.0064   0.000352842   0.000820158      0.9885      0.9885       85    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##  lp__       -20.01
##  m[1]         0.32
##  m[2]         0.44
##  cov[1,1]     1.04
##  cov[2,1]     0.56
##  cov[1,2]     0.56
##  cov[2,2]     1.26
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc,ch=F)
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##  lp__     -21.26 -20.89 1.83 1.55 -24.91 -19.06 1.00      674      854
##  m[1]       0.33   0.33 0.29 0.28  -0.14   0.80 1.00     1242     1225
##  m[2]       0.45   0.45 0.32 0.31  -0.05   0.96 1.00     1234     1121
##  cov[1,1]   1.61   1.43 0.70 0.49   0.86   2.98 1.00     1636     1341
##  cov[2,1]   0.90   0.79 0.64 0.44   0.18   1.98 1.00     1118      721
##  cov[1,2]   0.90   0.79 0.64 0.44   0.18   1.98 1.00     1118      721
##  cov[2,2]   1.98   1.78 0.96 0.66   1.02   3.62 1.00     1243      854

distributions

normal dist.

distribution of mean, it's from central limit theorem
y~N(m,s), y(-Inf,Inf), E[y]=m, V[y]=s^2

ex1-0.stan

data{
  int N;
  vector[N] y;
}
parameters{
  real m;
  real<lower=0> s;
}
model{
  y~normal(m,s);
}
generated quantities{
  real y1;
  y1=normal_rng(m,s);
}
n=20
y=rnorm(n,2,1)
summary(y)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   -0.22    0.99    1.65    1.77    2.30    3.79
par(mfrow=c(1,2))
hist(y,main='y')
density(y) |> plot(main='y')

data=list(N=n,y=y)

mdl=cmdstan_model('./ex1-0.stan')

mle=mdl$optimize(data=data)
## Initial log joint probability = -32.4913 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       10      -10.1074   0.000263813   0.000236765           1           1       16    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##      lp__   -10.11
##      m        1.77
##      s        1.01
##      y1      -0.25
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc)
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##      lp__ -11.08 -10.77 0.98 0.69 -13.06 -10.14 1.00      814     1047
##      m      1.78   1.78 0.25 0.24   1.38   2.18 1.00     1020     1006
##      s      1.10   1.08 0.19 0.17   0.85   1.45 1.00     1897     1416
##      y1     1.79   1.81 1.15 1.10  -0.12   3.71 1.00     1944     1884

Bernoulli dist.

The event with probablity p occur y=1, not occur y=0 
y~Ber(p), y{0,1}, E[y]=p, V[y]=p(1-p)

ex1-1.stan

data{
  int N;
  array[N] int y;
}
parameters{
  real<lower=0,upper=1> p;
}
model{
  y~bernoulli(p);
}
generated quantities{
  real s;
  s=(p*(1-p))^.5;
}
n=10
y=sample(0:1,n,replace=T)
summary(y)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##     0.0     0.0     0.5     0.5     1.0     1.0
par(mfrow=c(1,2))
hist(y,main='y')
density(y) |> plot(main='y')

data=list(N=n,y=y)

mdl=cmdstan_model('./ex1-1.stan')

mle=mdl$optimize(data=data)
## Initial log joint probability = -8.83422 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##        5      -6.93147    0.00017296    2.8168e-10           1           1        8    
## Optimization terminated normally:  
##   Convergence detected: gradient norm is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##      lp__    -6.93
##      p        0.50
##      s        0.50
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc)
##  variable  mean median   sd  mad     q5   q95 rhat ess_bulk ess_tail
##      lp__ -8.83  -8.55 0.72 0.32 -10.24 -8.32 1.00     1247     1431
##      p     0.51   0.51 0.14 0.14   0.27  0.73 1.01      817     1193
##      s     0.48   0.49 0.03 0.01   0.43  0.50 1.00     1248     1431

geometric dist.

The event with probablity p occur after y trials(failure) 
y~Ge(p), y[0,Inf), P(Y=y)=p(1−p)^y
E[y]=1/p, V[y]=(1-p)/p^2

ex1-2.stan

data{
  int N;
  array[N] int y;
}
parameters{
  real<lower=0,upper=1> p;
}
model{
  for (n in 1:N) {
    target+=log(p)+y[n]*log(1-p);
  }
}
generated quantities{
  real m;
  m=1/p;
  real s;
  s=((1-p)/p^2)^.5;
}
n=20
y=rgeom(n,0.3)
summary(y)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.00    1.00    2.00    3.05    4.00    9.00
par(mfrow=c(1,2))
hist(y,main='y')
density(y) |> plot(main='y')

data=list(N=n,y=y)

mdl=cmdstan_model('./ex1-2.stan')

mle=mdl$optimize(data=data)
## Initial log joint probability = -72.8551 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##        5      -45.2724    0.00138606   0.000172009           1           1        7    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##      lp__   -45.27
##      p        0.25
##      m        4.05
##      s        3.51
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc)
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##      lp__ -47.42 -47.17 0.65 0.30 -48.74 -46.95 1.00      903     1451
##      p      0.25   0.25 0.05 0.05   0.18   0.33 1.00      747      952
##      m      4.08   3.94 0.79 0.71   3.03   5.56 1.00      747      952
##      s      3.55   3.40 0.80 0.72   2.48   5.03 1.00      747      952

negative binomial dist.

The event with probablity p occur a times after y-a trials 
y~NB(a,p), y[a,Inf), P(Y=y)=p^a*(1−p)^(y-a)
E[y]=a*(1-p)/p, V[y]=a*(1-p)/p^2

ex1-3.stan

data{
  int N;
  array[N] int y;
}
parameters{
  real<lower=0,upper=1> p;
  real<lower=1> a;
}
model{
  y~neg_binomial(a,p);
}
generated quantities{
  real m;
  m=a*(1-p)/p;
  real s;
  s=(a*(1-p)/p^2)^.5;
  int y1;
  y1=neg_binomial_rng(a,p);
}
n=50
a=3
y=rnbinom(n,a,0.3)
summary(y)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##     0.0     3.0     5.0     6.1     8.0    22.0
par(mfrow=c(1,2))
hist(y,main='y')
density(y) |> plot(main='y')

data=list(N=n,y=y)

mdl=cmdstan_model('./ex1-3.stan')

mle=mdl$optimize(data=data)
## Initial log joint probability = -143.486 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       13      -138.138   0.000288552    0.00010991      0.9923      0.9923       20    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##      lp__  -138.14
##      p        0.43
##      a        2.65
##      m        3.45
##      s        2.82
##      y1       4.00
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc)
##  variable    mean  median   sd  mad      q5     q95 rhat ess_bulk ess_tail
##      lp__ -140.10 -139.74 1.17 0.84 -142.58 -138.98 1.01      530      705
##      p       0.52    0.50 0.16 0.15    0.29    0.84 1.03      344      436
##      a       3.14    3.02 0.95 0.90    1.83    4.92 1.03      364      410
##      m       2.96    2.99 1.11 1.04    0.92    4.70 1.02      430      447
##      s       2.50    2.45 0.87 0.80    1.04    3.97 1.03      373      447
##      y1      6.14    5.00 4.40 4.45    1.00   14.00 1.00     1826     1870

y1=mcmc$draws('y1')
par(mfrow=c(1,2))
hist(y1,main='y1')
density(y1) |> plot(main='y1') 

binomial dist.

The event with probablity p occur y times on n trials
y~B(n,p), y{0,1...n}, E[y]=np, V[y]=np(1-p)

ex2-3.stan

data{
  int N;
  int n;
  array[N] int y;
}
parameters{
  real<lower=0,upper=0> p;
}
model{
  y~binomial(n,p);
}
generated quantities{
  real m;
  m=n*p;
  real s;
  s=(n*p*(1-p))^.5;
  int y1;
  y1=binomial_rng(n,p);
}
n=20
n0=10
y=rbinom(n,n0,0.3)
summary(y)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    1.00    2.75    4.00    3.50    4.00    5.00
par(mfrow=c(1,2))
hist(y,main='y')
density(y) |> plot(main='y')

data=list(N=n,n=n0,y=y)

mdl=cmdstan_model('./ex2-3.stan')

mle=mdl$optimize(data=data)
## Initial log joint probability = -130.246 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##        3      -129.489    0.00238177   0.000593422      0.9786      0.9786        5    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##      lp__  -129.49
##      p        0.35
##      m        3.50
##      s        1.51
##      y1       5.00
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc)
##  variable    mean  median   sd  mad      q5     q95 rhat ess_bulk ess_tail
##      lp__ -131.47 -131.19 0.73 0.30 -132.88 -130.97 1.00     1055     1371
##      p       0.35    0.35 0.03 0.03    0.30    0.41 1.00      864     1086
##      m       3.52    3.52 0.34 0.33    2.96    4.07 1.00      864     1086
##      s       1.51    1.51 0.03 0.03    1.44    1.55 1.00      864     1086
##      y1      3.48    3.00 1.53 1.48    1.00    6.00 1.00     2011     2035

y1=mcmc$draws('y1')
par(mfrow=c(1,2))
hist(y1,main='y1')
density(y1) |> plot(main='y1') 

ex2-31.stan

trial n is varied from each sample

data{
  int N;
  array[N] int n;
  array[N] int y;
}
parameters{
  real<lower=0> p;
}
model{
  y~binomial(n,p);
}
n=20
n0=floor(runif(n,1,10))
y=rbinom(n,n0,0.3)

data=list(N=n,n=n0,y=y)

mdl=cmdstan_model('./ex2-31.stan')

mle=mdl$optimize(data=data)
## Initial log joint probability = -63.8628 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##        5      -57.6889   0.000555619   8.10574e-06           1           1        7    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##      lp__   -57.69
##      p        0.28
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc)
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##      lp__ -59.83 -59.55 0.75 0.34 -61.34 -59.30 1.00     1102     1428
##      p      0.28   0.28 0.05 0.05   0.21   0.36 1.00      711     1031

Poisson dist.

The event occur l times in unit range, the event occur y times. 
y~Po(l), y{0,1...Inf}, E[y]=l, V[y]=l

ex1-4.stan

data{
  int N;
  array[N] int y;
}
parameters{
  real<lower=0> l;
}
model{
  y~poisson(l);
}
generated quantities{
  real s;
  s=l^.5;
  int y1;
  y1=poisson_rng(l);
}
n=20
y=rpois(n,3)
summary(y)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##     1.0     2.0     3.0     3.1     4.0     6.0
par(mfrow=c(1,2))
hist(y,main='y')
density(y) |> plot(main='y')

data=list(N=n,y=y)

mdl=cmdstan_model('./ex1-4.stan')

mle=mdl$optimize(data=data)
## Initial log joint probability = -84.8157 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##        4       8.14693    1.7062e-05   2.66629e-05      0.9754      0.9754       10    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##      lp__     8.15
##      l        3.10
##      s        1.76
##      y1       2.00
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc)
##  variable mean median   sd  mad   q5  q95 rhat ess_bulk ess_tail
##      lp__ 8.80   9.05 0.69 0.33 7.51 9.28 1.00     1075     1174
##      l    3.15   3.14 0.39 0.40 2.54 3.81 1.00      835      947
##      s    1.77   1.77 0.11 0.11 1.59 1.95 1.00      835      947
##      y1   3.08   3.00 1.78 1.48 1.00 6.00 1.00     1925     1895

y1=mcmc$draws('y1')
par(mfrow=c(1,2))
hist(y1,main='y1')
density(y1) |> plot(main='y1') 

exponential dist.

The event occur l times in unit range, the event take time y to occur
y~ex(l), y>0, E[y]=1/l, V[y]=1/l^2

ex1-5.stan

data{
  int N;
  vector[N] y;
}
parameters{
  real<lower=0> l;
}
model{
  y~exponential(l);
}
generated quantities{
  real m;
  m=1/l;
  real s;
  s=1/l;
  real y1;
  y1=exponential_rng(l);
}
n=20
y=rexp(n,2)
summary(y)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.014   0.094   0.349   0.476   0.829   1.203
par(mfrow=c(1,2))
hist(y,main='y')
density(y) |> plot(main='y')

data=list(N=n,y=y)

mdl=cmdstan_model('./ex1-5.stan')

mle=mdl$optimize(data=data)
## Initial log joint probability = -20.9497 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##        2       -5.1408    0.00473542   1.33923e-09      0.5018           1        6    
## Optimization terminated normally:  
##   Convergence detected: gradient norm is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##      lp__    -5.14
##      l        2.10
##      m        0.48
##      s        0.48
##      y1       0.41
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc)
##  variable  mean median   sd  mad    q5   q95 rhat ess_bulk ess_tail
##      lp__ -4.89  -4.59 0.70 0.31 -6.38 -4.38 1.00     1089     1408
##      l     2.20   2.16 0.48 0.47  1.45  3.05 1.00      720     1189
##      m     0.48   0.46 0.11 0.10  0.33  0.69 1.00      720     1189
##      s     0.48   0.46 0.11 0.10  0.33  0.69 1.00      720     1189
##      y1    0.47   0.30 0.53 0.32  0.02  1.48 1.00     1860     1929

y1=mcmc$draws('y1')
par(mfrow=c(1,2))
hist(y1,main='y1')
density(y1) |> plot(main='y1') 

gamma dist.

The event occur l times in unit range, the event take time y to occur a times
y~Ga(a,l), y>0, E[y]=a/l, V[y]=a/l^2

ex1-6.stan

data{
  int N;
  vector[N] y;
}
parameters{
  real<lower=0> a;
  real<lower=0> l;
}
model{
  y~gamma(a,l);
}
n=20
y=rgamma(n,2,1)
summary(y)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.47    1.30    1.86    2.00    2.40    4.42
par(mfrow=c(1,2))
hist(y,main='y')
density(y) |> plot(main='y')

data=list(N=n,y=y)

mdl=cmdstan_model('./ex1-6.stan')

mle=mdl$optimize(data=data)
## Initial log joint probability = -37.5795 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##        7      -28.0203   0.000218728    0.00136526           1           1        9    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##      lp__   -28.02
##      a        3.33
##      l        1.67
##      m        2.00
##      s        1.09
##      y1       1.05
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc)
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##      lp__ -27.12 -26.84 0.91 0.73 -28.92 -26.19 1.00      627     1048
##      a      3.73   3.62 1.05 1.08   2.21   5.65 1.01      440      709
##      l      1.90   1.84 0.57 0.58   1.07   2.91 1.00      443      701
##      m      2.00   1.97 0.24 0.22   1.65   2.43 1.00     1811     1184
##      s      1.06   1.04 0.20 0.19   0.79   1.43 1.00      525      784
##      y1     2.04   1.85 1.12 0.98   0.65   4.08 1.00     1979     1932

y1=mcmc$draws('y1')
par(mfrow=c(1,2))
hist(y1,main='y1')
density(y1) |> plot(main='y1') 

log normal dist.

logalithm of variable follows normal distribution
log y~N(m,s), y>0
E[y]=exp(m+s^2)
V[y]=exp(2*m+s^2)*(exp(s^2)-1)

ex1-7.stan

data{
  int N;
  vector[N]<lower=0> y;
}
parameters{
  real m0;
  real<lower=0> s0;
}
model{
  y~lognormal(m0,s0);
}
n=30
y=rlnorm(n,0.5,1)
summary(y)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.25    0.97    1.66    3.75    4.38   23.99
par(mfrow=c(1,2))
hist(y,main='y')
density(y) |> plot(main='y')

data=list(N=n,y=y)

mdl=cmdstan_model('./ex1-7.stan')

mle=mdl$optimize(data=data)
## Initial log joint probability = -75.9714 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##        7      -44.0063   8.10734e-05   6.20292e-05           1           1       12    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##      lp__   -44.01
##      m0       0.74
##      s0       1.05
##      m        6.32
##      s        5.16
##      y1       6.81
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc)
##  variable   mean median    sd  mad     q5    q95 rhat ess_bulk ess_tail
##      lp__ -44.97 -44.67  1.01 0.72 -47.07 -44.00 1.00      831      992
##      m0     0.74   0.74  0.21 0.20   0.39   1.08 1.00      818      882
##      s0     1.11   1.10  0.15 0.14   0.90   1.39 1.00     1891     1522
##      m      8.09   7.06  4.22 2.47   4.29  15.49 1.00     1376     1254
##      s      6.92   5.87  4.14 2.38   3.26  14.16 1.00     1491     1491
##      y1     4.22   2.07 10.10 2.04   0.34  13.78 1.00     2001     2055

y1=mcmc$draws('y1')
par(mfrow=c(1,2))
hist(y1,main='y1')
density(y1) |> plot(main='y1') 

beta dist.

use as prior of binomial dist parameter p
event with probabilty p0 occur a-1 times, do not occur b-1 times in a+b-2 trials
p~Be(a,b), p[0,1], p=p0^(a-1)*(1-p0)^(b-1)
E[p]=a/(a+b), V[p]=ab/(a+b)^2/(a+b+1)

ex1-8.stan

data{
  int N;
  vector[N] p;
}
parameters{
  real<lower=0> a;
  real<lower=0> b;
}
model{
  p~beta(a,b);
}
n=20
p=rbeta(n,1,2)
summary(p)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.025   0.149   0.250   0.316   0.446   0.721
par(mfrow=c(1,2))
hist(p,main='p')
density(p) |> plot(main='p')

data=list(N=n,p=p)

mdl=cmdstan_model('./ex1-8.stan')

mle=mdl$optimize(data=data)
## Initial log joint probability = -60.7305 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       11        6.3322   0.000267379   0.000143815      0.9735      0.9735       13    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##      lp__     6.33
##      a        1.48
##      b        3.24
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc)
##  variable mean median   sd  mad   q5  q95 rhat ess_bulk ess_tail
##      lp__ 6.96   7.29 1.05 0.79 4.81 8.00 1.00      848      981
##      a    1.67   1.61 0.47 0.46 0.97 2.51 1.01      528      654
##      b    3.70   3.61 1.13 1.14 1.98 5.63 1.01      542      793

Dirichlet dist.

use as prior of categorical dist parameter p
p~dir(p0), p[0,1]

ex1-9.stan

data {
  int N;
  int K;
  matrix[N, K] p;
}
parameters {
  simplex[K] p0;
}
model {
  for(i in 1:N){
     p[i]~dirichlet(p0);
  }
}
library(gtools)
n=20
p=rdirichlet(n,c(0.5,0.3,0.2))
summary(p)
##        V1              V2              V3       
##  Min.   :0.001   Min.   :0.000   Min.   :0.000  
##  1st Qu.:0.183   1st Qu.:0.012   1st Qu.:0.003  
##  Median :0.842   Median :0.098   Median :0.034  
##  Mean   :0.615   Mean   :0.299   Mean   :0.085  
##  3rd Qu.:0.927   3rd Qu.:0.617   3rd Qu.:0.085  
##  Max.   :0.998   Max.   :0.919   Max.   :0.724
boxplot(p)

data=list(N=n,K=ncol(p),p=p)

mdl=cmdstan_model('./ex1-9.stan')

mle=mdl$optimize(data=data)
## Initial log joint probability = 76.6785 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##        5       81.6749   0.000231235    6.7112e-05           1           1        8    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##     lp__     81.67
##     p0[1]     0.55
##     p0[2]     0.26
##     p0[3]     0.18
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.3 seconds.
seeMCMC(mcmc)
##  variable  mean median   sd  mad    q5   q95 rhat ess_bulk ess_tail
##     lp__  77.03  77.32 1.00 0.75 75.05 78.01 1.00     1011     1268
##     p0[1]  0.54   0.54 0.06 0.06  0.45  0.64 1.00     1794     1473
##     p0[2]  0.27   0.26 0.05 0.05  0.19  0.35 1.00     1795     1223
##     p0[3]  0.19   0.19 0.04 0.04  0.13  0.25 1.00     1844     1453